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The real time evolution of two pieces of quantum insulators, initially at different temperatures, is 
studied when they are glued together. Specifically, each subsystem is taken as a Bose-Hubbard model 
in a Mott insulator state. The process of temperature equilibration via heat transfer is simulated 
in real time using the Minimally Entangled Typical Thermal States algorithm. The analytic theory 
based on quasiparticles transport is also given. 
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I. INTRODUCTION 

One of the open problems of statistical mechanics is 
the description of isolated quantum many-body systems 
far from equilibrium. Simulations of real time dynam¬ 
ics are quite difficult and requiring large computer re¬ 
sources even for the simplest nontrivial one-dimensional 
systems with short range interactions, although a signif¬ 
icant progress has been obtained in the last ten or so 
years when modern versions of Density Matrix Renor¬ 
malization Group (DMRG) [1] approaches using Matrix 
Product States (MPS) have been developed [2, 3] - for 
recent reviews see [4, 5]. It has been understood that the 
growth of entanglement in the system is a main obsta¬ 
cle in studying long-time dynamics even for initial states 
being close to the ground state, for which the entangle¬ 
ment is typically quite low [5] . This growth may be even 
more restrictive for abrupt quenches and study of excited 
states [6]. 

By comparison, the study of finite temperature dynam¬ 
ics is much less understood. Some results for large sys¬ 
tems were obtained assuming adiabatic evolution (with 
entropy conservation [7, 8]). One possible “quasi-exact” 
approach involves MPS of a system enlarged by an aux¬ 
iliary environment. By tracing out the auxiliary degrees 
of freedom, one obtains the density matrix of the original 
system [9, 10]. This approach has been applied for dy¬ 
namics of thermal systems [11-13]. Interestingly, there 
is an important freedom in the approach, namely that 
the same density matrix of the system may be obtained 
for different dynamics of the environment. One may use 
this freedom in the wavefunction describing the system 
and the environment, in order to minimize the growth of 
system-environment entanglement during the temporal 
evolution. That in turn allows for reaching longer times 
of evolution [14, 15]. This approach has been used to 
calculate particle currents and Drude weight as well as 
spin correlations in XXZ spin chains. 

Another promising approach is based on the combina¬ 
tion of time-dependent DMRG in the Heisenberg picture 
- that allows to obtain the real time evolution of the 


operator [16] - with the density matrix obtained using 
imaginary time propagation. This allows for an evalua¬ 
tion of expectation values on a grid of temperature/time 
points simultaneously [17]. 

We shall use in the present paper yet another ap¬ 
proach called Minimally Entangled Typical Thermal 
State (METTS) approach [18-20]. It has been argued 
that it can simulate finite temperature quantum systems 
with a computational cost comparable to ground state 
DMRG [19]. As its application to real time dynamics 
requires propagation of excited states for which the en¬ 
tanglement growth may be a serious obstacle [21], this 
method may be costly to implement. A very recent study 
[22] has shown that METTS may be quite efficient for 
gaped systems at low temperatures. We shall use it for 
studying the dynamics of gaped Mott insulators far from 
equilibrium. 

We will mainly consider a system composed of two one¬ 
dimensional Mott insulators at different temperatures 
that are glued together at a given instant of time. Such 
a system can, in principle, be realized in the laboratory. 
A single insulator is routinely observed by placing ultra¬ 
cold atoms in a one-dimensional deep optical lattice with 
the transverse directions being frozen by a tight laser 
confinement [23]. The system can be split by an addi¬ 
tional laser potential into two separate parts. One may 
imagine heating both parts differently and then bringing 
them into contact by a rapid switch-off of the separating 
laser. 

For such a model system, we concentrate on the heat 
transfer assuming little direct particle current. In Section 
II, we discuss briefly our METTS implementation for the 
Bose-Hubbard Hamiltonian and show exemplary results. 
Section III brings a simple analytical theory for the ob¬ 
servables assuming that the heat transport is dominated 
by quasiparticle motions, described (in the low tunneling 
limit) by the Bogoliubov approach. Both approaches are 
compared in Section IV where we also discuss an alterna¬ 
tive approach in which the two subsystems are smoothly 
glued. We conclude in Section V. 
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II. THE SYSTEM STUDIED USING METTS 


B. System preparation 


A. METTS algorithm 


The METTS algorithm as proposed by White [18] sim¬ 
ulates a thermal canonical ensemble with temperature 
T = 1//3 (from now on we shall assume the convenient 
units with the Boltzmann constant kp = 1, we assume 
also h = 1). Since the method is described in detail 
elsewhere [19], we provide essentials only. The METTS 
approach works by alternately applying the imaginary 
time evolution operator exp(— fiH/2) and a projection 
measurement onto a given basis set. This defines a ran¬ 
dom walk which samples the Hilbert space and creates 
a properly weighted ensemble M which enables to esti¬ 
mate thermal averages of operators in the canonical en¬ 
semble: (O) = Tr[exp (—/3H)0\ as (O) = Y 

peM 

The METTS ensemble allows also for real time evo¬ 
lution of the ensemble ( 0(t )) = Y ( / 0|O(t)|'0) = 

Y ($(t)\0\^(t)) . Both the imaginary time and the real 

time evolutions are performed efficiently with the time- 
evolving block decimation (TEBD) algorithm [2], essen¬ 
tially equivalent to a time-dependent DMRG approach. 

The Hamiltonian of the one-dimensional Bose- 
Hubbard (BH) model reads 


i i= 1 K 




jj > 

Oi4 +1 + h.c + —riiiifii - 1 ) 


(1) 

with a,i (aj) being the standard boson annihilation (cre¬ 
ation) operator at site i = 1,...,L, rii = aja*, R = J is 
the amplitude of hopping (tunneling) between i and i +1 
site; Jl = 0 for open boundary conditions (OBC) consid¬ 
ered below while U denotes the interaction strength. We 
consider the insulating regime V J (we shall assume 
typically that J/U = 0.05 deep in the Mott regime). In 
the METTS approach, the projection measurement can 
be performed on any basis set. It is here convenient to 
use a Fock basis on each lattice site, meaning that we 
measure the (random) number of particle on each site, 
an easy task for a MPS state new [24]. For low temper¬ 
atures considered in this work (we take a typical /3U of 
the order of 5), successive evolution steps between mea¬ 
surements are correlated. To avoid that, we have started 
the METTS algorithm from the ground state neglecting 
the first 500 iterations of the METTS sampling proce¬ 
dure. For the final ensemble, we have taken every 200th 
METTS vector. Such a sampling is quite time consum¬ 
ing. The alternative would be to change the basis in 
which the measurements are performed. That, however, 
would significantly slow down the time evolution using 
TEBD, because it would break the total particle number 
conservation (which would of course be restored after sta¬ 
tistical averaging) [2]. 


Our aim is to study the transport in a non-equilibrium 
system of two insulators at different temperatures 
brought into contact at a given time. Two versions of 
the procedure will be considered: 

• Two uniform systems with different temperatures 
are merged by building the tensor product of their 
density matrices. 

• Preparation of a temperature-inhomogeneous sys¬ 
tem with a smooth temperature gradient across the 
sample, using the canonical thermal distribution of 
an auxiliary Hamiltonian. 


1. Simple tensor product approach 

To prepare the initial situation, we have considered two 
lattices of the same length M = 35 sites that are linked 
together forming a longer lattice of length L = 2M. The 
Hamiltonian is H = Hp + Hr + Vlr where both Hi 
terms take the form (1) and Vlr is the hopping term 
connecting the two lattices. Let us define a nonstan¬ 
dard but useful convention that the two systems are 
linked “symmetrically” at i = 0. Thus the indices of 
the right hand side lattice of M sites take half integer 
values i = 1/2, 3/2,..., M/2 — 1/2 while those of the left 
hand side lattice are i = —M/2 + 1/2,... — 3/2,—1/2. 
Then the coupling term between two subsystems reads 
— J(cl_i/ 2 Cl \/2 + h.c.). We shall assume this “symmetric” 
notation from now on. 

The density matrix that describes the full system 
would be just the product of constituents’ density ma¬ 
trices: p = pL 0 each prepared separately un¬ 
der open boundary conditions at different temperatures, 
i.e., pi = exp(— PiHi) with i = L,R. We compute 
two appropriate METTS Ml and Mr that represent 
the density matrices pL and pR. The METTS ensem¬ 
ble that represents the density matrix pp 0 Pr is just 
Ml 0 Mr := {Vt 0 ^>r\A G Mi,i = L,Hj. This is 
not efficient computationally. For example when estimat¬ 
ing (Oi) for i e {L,R} the actual average is as follows 
(in this case for L): (^L|Oi|Vh,)(fMV’ij) = {fpL\Oi\ip L ). 
The sum over METTS vectors from Ml 0 Mr will 
contain repetitions of each average. It will occur as 
many times as many vectors are in Mr. On the other 
hand, if the two METTS ensembles contain exactly the 
same number of vectors V, then the METTS ensem¬ 
ble Ml®Mr yields completely equivalent estimates as: 
Ml 0 Mr := {^p 0 ^ l R \i = 1,..., U}, where ^ l L/R de¬ 
notes the i-th METTS vector of Ml/r. 

In typical applications to the BH model, we have found 
that the size of Ml 0 Mr is several millions of METTS 
which is computationally prohibitive, while Ml 0 Mr 
contains only the same number of vectors as each of Ml 
and Mr. We have therefore used the second approach. 
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FIG. 1. Dynamical evolution of the variance of the particle 
number on each site. The initial state is chosen with an in¬ 
homogeneous temperature profile, in order to induce an heat 
current. The black curves are obtained for a tensor product 
initial state, where the left and right parts of the system are 
both at thermal equilibrium, but at different temperatures 
PlU — 6, /3rU m 8. The red curves correspond to "smoothly" 
glued initial conditions as described in the main text. For 
tU = 2, the cusp in the black curve is a remnant of boundary 
effects occurring when tunneling between the two subsystems 
is abruptly turned on at t — 0. The short range fluctuations 
on the curves are of statistical origin, an unavoidable effect 
of the METTS approach. They vanish when the number of 
METTS states tends to infinity. 


The Hamiltonian governing the dynamics of the whole 
lattice couples the “left” and “right” parts of the system 
through the additional tunneling term, mentioned above, 
— Ja_ 1 / 2 a\/ 2 + h.c.. When using the tensor product of 
two density matrices as an initial state, the gluing region 
is subject to violent short-time transient phenomena as 
shown in Fig. 1: see for example the increased variance 
just on the left of the gluing point at tU = 2. They are 
the remnants of the OBC for the two constituent lattices. 
This region is also the region where temperature gradient 
should be expected. Note that in particular this merging 
method implies that the notion of "temperature” in the 
middle of the system is questionable. On the other hand 
this approach seems the simplest one. It will be referred 
to as the tensor product approach. 

The data shown in Fig. 1 - as well as in all results 
of quasi-exact METTS simulations in the sequel of this 
paper - display site-to-site fluctuations. This is an in¬ 
trinsic drawback due to the statistical description of the 
system state in the METTS method; it also affects all 
statistical methods a la Monte-Carlo. These short range 
fluctuations decrease when the number of METTS states 
increases. Note also that they decay in the course of the 
temporal evolution, as clearly seen in Fig. 1. They are 
responsible for the noisy character of some figures shown 
below, but they do not affect any of the conclusion of this 
paper. 


2. Smooth Gluing 

Another approach is possible which makes the tran¬ 
sition region between the two subsystems more subtle. 
We prepare the initial state in just a single step: the 
inhomogeneous system containing the “left” and “right” 
parts with different temperatures is prepared in a sin¬ 
gle canonical ensemble simulation at thermal equilibrium. 
To achieve that, we consider an auxiliary Hamiltonian 
T-L = JA with the scaling function / smoothly in¬ 

terpolating between ^ on the left side and /3r on the 
right side. For example, we have used: 

r/.x Pr-Pl, , ( ., x . Pr + Pl 
/(*) = -g - tanh (*/«) +-g - ( 2 ) 

where 1 s <C M denotes the size of the interface. 
We then construct the canonical density matrix for the 
auxiliary Hamiltonian at (3 = 1. Then, provided the 
correlation length in the system is much shorter than the 
interface size s it is expected that the reduced density ma¬ 
trix of the left part pL = Tr^exp (—T-L) ~ exp(— PlHl), 
while that for the right part reads Trz,exp(— Ti) « 
exp(— (3rHr). This agreement requires also achieving a 
thermodynamic limit in terms of system size. 


C. Observables and thermometry 

The typical observables measured in the created en¬ 
semble or during the subsequent temporal evolution is 
the average number of particles per site (n$) as well as 
its variance Var(nQ = (n?) — (nQ 2 . The variance at 
finite temperature contains both the quantum and the 
thermal fluctuations and can be used as a measure of the 
temperature. The advantage of this measure is that it 
is local, and thus well suited for the situation of interest 
where there are a space-dependent temperature profile 
and heat currents. It is however a priori not obvious 
that the local variance is a faithful measure of the local 
temperature. In order to test this assumption (similar in 
spirit to the well-known local density approximation used 
for optical lattices exposed to a slowly spatially varying 
trapping potential, see e.g. [25]) we consider smoothly 
glued samples - see section IIB 2 with M = 35 (at unit 
mean density: TO bosons on TO sites), J/U = 0.05 and 
P L U = 8 , p R U = 12 . 

The Quantum Monte Carlo (QMC) approach as im¬ 
plemented in ALPS [26] allows to find the canonical en¬ 
semble density of T-L : p oc exp (—T-L) and consequently to 
compute the the variance Var(nQ for a temperature in¬ 
homogeneous system. The function /(i), eq. (2), is the 
local (inverse) temperature. In Fig. 2, we show the lo¬ 
cal variance Var(nQ vs. the local inverse temperature. 
Except for small differences related to finite size effects, 
all the data collapse on the same curve, regardless of the 
sharpness of the transition between the hot and the cold 
parts of the system. 
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FIG. 2. (color online) Local thermometer based on the parti¬ 
cle number variance, (a): local variance on the various sites 
for two smoothly glued Mott insulators, with inverse temper¬ 
atures [3lU = 8 (left part) and /3rU = 12 (right part). The 
temperature profile is given by Eq. (2), with a (smooth) in¬ 
terface size s = 3 (black solid line), s = 7 (red, dashed line) 
and s = 10 (blue, dot-dashed line), (b): local variance vs. 
the (inverse) local temperature, Eq. (2). Except for small 
differences, the three curves collapse, proving that the local 
variance - with a good approximation - depends only on the 
local temperature, regardless of the sharpness of the temper¬ 
ature profile. This proves that the local variance can be used 
as a reliable thermometer in a inhomogeneous system. 


D. Real time evolution 

Real time evolution of each member of a METTS en¬ 
semble is performed using a home-made implementation 
of the TEBD algorithm, taking advantage of total par¬ 
ticle number conservation [27]. The ensemble consists 
typically of about 12000 METTS. By restricting to MPS 
of maximum bond dimension x = 150, we have been able 
to perform numerical evolution up to time tU = 120, 
(corresponding to tJ = 6) keeping discarded weights at 
the level of 10 -4 at most. We have used a fourth order 
Trotter decomposition in time to control the time dis¬ 
cretization error [28]. As we are able, see sections IIB1 
and IIB 2, to prepare an initial state far from the ther¬ 
mal equilibrium, and we can measure both its density 
and temperature profiles as a function of increasing time, 
we can directly monitor heat and mass transport in the 
system. We can moreover study quantitatively transport 
properties by calculating the energy and particle currents 


in the system. We follow to some extent the approach 
presented in [15] for spin systems. 

Let us rewrite first the Bose-Hubbard system in a sym¬ 
metrized form as H = JT hi with 

h% — ~ T CLi-\-\CL^ T tti&i_\ T CLi—lCL^j 

+ ^ni(ni - 1). (3) 

The currents flowing in and out of site i may be defined 
via continuity equations. For example the energy current 
may be obtained from 

d t hi = i[H,hi\ 

i\hi — 2 •> h"i\ + i\hi —l^ h^ H - i[hi- |_i, hi] x\hi- 1_2 5 hi] 

= j E (i ~ 1/2) — j E (i + 1/2), (4) 

where j E (i — 1/2) is the current coming into site i from 
site i — 1. Since the site index i is in our case a half-integer 
number, the current index is an integer. In particular 
^(0) corresponds to a current passing through the center 
where the two subsystems are glued together. A little 
algebra shows that 

j (i 1/2) i\hi— 2 , hi] T i\hi— i, hi] T i[hi— i, i], (5) 

with the first and last terms being the consequence of 
our symmetrized form of hi in (3) which involves a site 
energy plus half of the links to neighbors on both sides. 

The Hamiltonian for short range interactions on a lat¬ 
tice may be split as a sum of site terms in several ways, 
leading to slightly different expressions for the currents. 
For example, [15] uses an asymmetric site Hamiltonian 
containing a site energy and the full link to the right. 
However, the different forms lead to very similar results 
whenever the system changes smoothly from site to site. 

For our Bose-Hubbard system, the energy current op¬ 
erator defined in Eq. (5) reads explicitly: 

j E (i - 1/2) = £ [£(* - 2, i) + K[i - 1, * + 1)] 

+ ^[V(M-1)-V(i-M)] (6) 

with JC(i,j) = i{cb\aj — a]a,) and M(i,j) = i(a\niCLj — 
h.c.) (with rii = a\a,i). 

Similarly we may define the mass (particle) current 
j n {i = 1/2) [29] using the continuity equation for dtrii 

d t rii = i[H,rii] 

= i[hi-i,rii] + i[hi,rii] + i[h i+1 ,rii\ 

U j n {i - 1/2)-j n {i +1/2) (7) 

yielding 

j n (i - 1/2) = JJC(i,i-l). (8) 

Having all the observables defined, let us take a closer 
look at exemplary numerical data. We consider two 
sharply connected Mott insulators (with M = 35 each) 
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FIG. 3. (color online) Dynamical evolution of the aver¬ 
age particle number (m) and average squared particle num¬ 
ber (vp ) after two Mott insulators at different temperatures 
PlU = 6,/3rU = 8 are sharply connected at time t = 0. 
Each subsystem has initially 35 particles on 35 sites. Note 
the spreading of the central perturbation observed both on 
the average density (m) (lower, black curves in all the pan¬ 
els) and on the average squared density (wp ) (higher, red lines 
in all the panels). 


at PlU = 6 and PrU = 8. Again the system contains 70 
particles in total. For short times, the variance for this 
run has already been shown in Fig. 1. Fig. 3 shows the 
average density, (n*), and the averaged squared density, 
(nf), for longer times. Firstly observe that the site to 
site fluctuations of both quantities are larger than on the 
variance {jnp) — (rii) 2 (subtraction reduces fluctuations), 
compare with the lowest panel in Fig. 1). Apart from sta¬ 
tistical fluctuations one may clearly observe the spread¬ 
ing of the perturbation from the center of the sample with 
time. The density shows a small but clear excess above 
one on the right hand side (the cold one), this excess 
moves to the right with an apparently constant velocity. 
This excess must be somewhere compensated by a lack 
of particles (the particle number is conserved): indeed a 
hole moves to the left (hot side). Similarly bumps and 
holes in (vp) spread out in both directions approximately 
linearly in time. 

This picture is confirmed by observing the energy cur¬ 
rents shown in Fig. 4. One expects that, when gluing 
two subsystems together, the effect appears first around 
the gluing point. Then the information travels within 
the system modifying various observables (including cur¬ 
rents). We observe this effect in Fig. 4. The current at 
the place of merging, ^(0) almost immediately picks up 
non zero values starting from the gluing time t = 0. A 
crucial observation is that, after some transient, the cur¬ 
rent ^(0) is constant at long time, although the spatial 
profile of (rti) and (up) becomes smoother and smoother. 



FIG. 4. (color online) Energy currents associated with the 
data of Fig. 3. The current initially starts flowing at the 
connecting point, i — 0, between the two samples at different 
temperatures. It then propagates on both sides, reaching sites 
i = ±5 around tU — 25, and site i — 10 later around tU = 50. 
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FIG. 5. (color online) Spatio-temporal development of the 
energy current for the data represented in Fig. 3. Observe 
the ballistic spreading of the perturbation, initially born at 
the connecting point, i — 0, between the two Mott insulators 
at different temperatures. See text for further discussion. 


This definitely rules out the possibility of diffusive en¬ 
ergy transport in the system. Indeed, in such a case, 
the current would have been proportional to the local 
temperature gradient, and would decay at long times. 
The persistence of a finite current at long times is a di¬ 
rect signature of ballistic transport. This is confirmed 
by looking at the current at other sites. They fluctu¬ 
ate around zero at short time, and rise only after some 
position-dependent delay. The spread of the information 
is approximately linear in time (ballistic) as visible in 
Fig. 4 and more apparent while inspecting the color plots 
in Fig. 5. The dashed lines give the predicted theoreti¬ 
cal limits - derived in the next Section - for the ballistic 
spreading of information. 
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III. BOGOLIUBOV THEORY 
A. Reduction of the Hilbert space 

For a moderate value of the tunneling amplitude J, the 
Bose-Hubbard ground state at integer filling is a Mott in¬ 
sulator [30], with a non-zero gap. In order to understand 
heat and mass transport at finite temperature, it is thus 
useful to study the properties of the ground state as well 
as those of the low energy excitations. In this section, we 
give a complete quantitative description of these quanti¬ 
ties, in the limit of small tunneling J U and relatively 
low temperature /3U 1. Note that we do not put any re¬ 

striction on the parameter /? J, so that this approach takes 
into account both the quantum fluctuations induced by 
J, the thermal fluctuations governed by /3, and mixed 
contributions. We show below that this simple theory is 
able to quantitatively reproduce the quasi-exact numer¬ 
ical results obtained for heat and mass transport in the 
fully out of equilibrium Bose-Hubbard model. 

In the limit of vanishing J, the ground state at unit fill¬ 
ing is simply a Fock state with exactly one particle per 
site. Elementary excitations are particle-hole excitations 
where a second particle is added on one site while leaving 
another site empty. This costs an energy U. If one drops 
the requirement of a fixed total number of particles - thus 
going from a canonical to a grand canonical ensemble - 
elementary excitations are independent particle or hole 
excitations. The energy cost of the particle excitation is 
U — /x, where /z is the chemical potential, and the cost of 
a hole excitation is fi. The system is thus gaped for any 
value of /i in the [0,1/] interval. At finite temperature, 
particles or holes can be created, with probabilities re¬ 
spectively and e _/3 U The average density will 

be maintained at 1 if /jl = U/2. Higher excitations (sev¬ 
eral particles and/or holes) are negligible if [3U <C 1. 

For a small but finite J, the ground state is no longer 
a Fock state, but it remains gaped. Using a perturbative 
approach in powers of J/U, it is possible [31] to obtain 
accurate estimates of all interesting quantities such as 
energy, average occupation number (jii) or (n?) on site i. 
We will here use a similar approach, restricting to lowest 
non-vanishing order in J/J7, but including also elemen¬ 
tary excitations. While particle/hole excitations on all 
sites are degenerate for J = 0, tunneling between sites 
lifts this degeneracy producing a band of quasiparticle 
and quasihole excitations. 

At lowest order in <//[/, it is enough to consider sin¬ 
gle particle and hole excitations, that is to restrict the 
Hilbert space on each site to occupation numbers rii = 
0,1,2. This problem has been addressed in differents 
works,see e.g. [32, 33]; for a recent review with an exten¬ 
sive list of references, see [34]. In our description below, 
following [35], we define a ’Vacuum” state as the Fock 
state with one particle on each site and introduce particle 
and hole creation operators on site i : c\ and d \. Their ac¬ 
tion is best illustrated by mapping the boson occupation 
number onto two quantum numbers (n c , rid) representing 


particle and hole occupation numbers: on each site, we 
have |2) = |(1,0)>, |1) = |(0,0)> and |0) = |(0,1)>. One 
then has to implement that there cannot be less than zero 
or more than two bosons on each site and that states with 
one particle and one hole on the same site are forbidden. 

This results in an effective Hamiltonian quadratic in 
particle/hole operators, see [35] for a detailed derivation: 


H eS = - J £ N c i + 4dj + Aadj + c\d]) 1 +U £ 


(ij) 


ct Cj 


(9) 

where (z, j) denotes a pair of neighboring sites. 

When the density of excitations is small, 
(cjci), (cjci) < 1, the Ci,cj,di,dl can be approxi¬ 

mately considered as standard bosonic annihilation and 
creation operators. 


B. Bogoliubov approach 


For a finite system with periodic boundary conditions 
(or an infinite system), this effective Hamiltonian can be 
diagonalized using the translational invariance, going to 
the momentum space [35]: 

Ck = 4j2c r e~ ikr , 4 = / E , re -- (10) 


where L is the number of sites of the system. The effec¬ 
tive Hamiltonian writes: 


Hq ff — 2 J 'y ^ cos k 2c k c k + d k d k \Z2(^c k d— k ~h v k d k ) 

+ UE k c{c k (11) 


which can be conveniently rewritten as: 


tfeff = £[4, 
k 


U — 4 J cos k 

—2\f2J cos k 

Ck 

2\/2J cos k 

2 J cos k 

A k . 


( 12 ) 


This Hamiltonian can be diagonalized using a Bogoliubov 
transformation: 


Ck 


= u k B k + v* k A ] _ k , d)_ k = v k B k + u* k A ] _ k (13) 


where (uk,Vk) is the eigenmode of the Bogoliubov-de 
Gennes equations: 


(qp) 


E 


with eigenenergy: 


^k 


U — 4 J cos k 

—2 \/2J cosk 

u k 

_ v k _ 


2\[2 J cos k 

2 J cos k 

_ Vk 


- J cos k + GJfc, 


14) 


(15) 


where 


VU 2 - 12 JU cos k + 4 J 2 cos 2 k 

Uk = -X- 


(16) 
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With the normalization \uk \ 2 — \vk\ 2 = 1, the operators 
Ak,Bk,Al, B k are standard bosonic annihilation and cre¬ 
ation operators. The other eigenmode (—Vk,Uk) is asso¬ 
ciated with eigenenergy: 

E j: qh ) = ~ T J cos k T ujk , (17) 

so that the Hamiltonian writes: 


#e ff = £ 


El qp) BtB k + E, 


(q h ) /it 


U 


A\A k + u> k + 3 J cos k - — 


(18) 

Eqs. (15-17) are special cases of general formulae avail¬ 
able for arbitrary integer mean density (see e.g. [33]) as 
reviewed by Krutisky [34], see his Eq. (325). 


C. Ground state 


The ground state of the system is the vac¬ 
uum of the Ak,Bk operators and its energy is 
(ujk + 3 J cos/c — U/ 2 ). In a finite system with peri¬ 
odic boundary conditions, the k values are discretized as 
integer multiples of 2i t/L. The limit of large system is 
obtained with the substitution: 



In the limit of small J we are interested in, one gets: 

U J 2 

uj k « — — 3 J cos k - 8 — cos 2 k (20) 

so that the ground state energy is approximately equal 
to —4 J 2 L/U, (—AJ 2 /U per site), in agreement with [31] 
(note that the first order in J cancels out). 

It is also possible to compute expectation values of 
simple operators in the ground state. The total number 
of particles is: 

N = Y, n i = E + c k - d \di) = L + J2 C W (21) 

i i k 

It is simply expressed as: 

N = L + J2(BlB k -A{A k ) (22) 

k 

so that B k (resp. A^ k ) appears as the creation operator 
for a quasiparticle (resp. quasihole) with momentum k. 
In the ground state, one trivially recovers (N) = L, so 
that there is on average exactly one particle per site. 

The variance of the number of particle per site is more 
interesting. Indeed: 


N 2 =Y J n i =Yj( 1+ C k “ d \ d i) 2 ’ 

i i 


which gives, in the limit of low density of excitations: 


N 2 =L + J2 (3 c\ Ci - d\di) = L + J2 (34 c fc - 44). 

i k 

(24) 

The physical interpretation of the +3 and -1 coefficients 
is rather transparent: c\ creates an additional particle 
at site i, thus n 2 = 4, i.e. 3 on top of the 1 particle of 
the background. Similarly, d\ produces rii = 0, that is 1 
below the background. 

In turn, the operator N 2 can be expressed as a function 
of the Ak and B operators: 


N 2 =L + E k 


(3M 2 


\ v k\ 2 )B k B k - (\u k \ 2 - 3|t)fc| 2 )A\A k 


+2u k v k A k B k + 2t i^v^A^B^ + 2 |v^p 


(25) 


In the limit of small J, the eigenmodes of Eq. (14) are: 


Uk = 1 , Vk = 2A/2COS k J/U , (26) 


which, in a more general form can be found also in 
[34], see Eq. (326) there. The expectation values on the 
ground state is: 

t 2 / o t2 \ 

(AT 2 ) = L + ^16cos 2 fc^ = L(l+^), (27) 

k ^ ' 


recovering the fact that the variance of the number of 
particle per site is 8 J 2 /U 2 [31]. 


D. Thermal excitations 


We now turn to thermal excitations in the system. One 
can use either the canonical or the grand canonical en¬ 
semble, all expectation values are expected to be equal in 
the two ensembles in the large L limit. It turns out that 
the grand canonical ensemble is more convenient for cal¬ 
culations. We thus have to consider configurations with a 
weight scaling like The operator H e q — fiN 

is readily obtained from Eqs. (18),(22): 


ffeff ~^N = [( 4 qp) - v)B\B k + ( 4 qh) + ^)A\A k 

+c dk + 3 J cos k — - 7 ^] . (28) 


The calculation of expectation values at the thermal equi¬ 
librium is thus straightforward. One obtains (Bk) = 


(• B l ) = ( A k ) = (-41) = 0 and: 


Tr 

'BlB k e-^ E * v) -^ B l 

B k ~ 

Tr 

e-/3(4 qp) -M)4 B *' 



(29) 


which, in the limit of low density of excitations, reduces 
to: 






( 23 ) 


( 30 ) 












Similarly: 


<4 A k ) 


p -/3 (4“ h) +^) 


(31) 


We thus obtain for the average number of particles: 


(TV) 




-/?(S< qp) -M) _ „-/3(E^ qn; +M) 


(qh)_ 


(32) 


Again, the interpretation in terms of density of quasi¬ 
particles and quasiholes is simple as: 

^ = 1 + 8 ^+ f +7T [3n^\k)-n^ h \k)} ( 40 ) 

L U z J_ 7r 2 tt 

The variance is finally obtained by combining this 
equation with Eq. (36): 


In the limit of small J/U, one can use the approximate 
expression (20) at first order and obtain: 


(N) = L + 


^ ^ /3(£/—4 J cos k—fi) 

k 


cos k ) 


(33) 

In the continuous limit of large U, the sum over k becomes 
an integral giving: 


4 = 1 + e-^-^XotW) - e“ /3A1 X 0 (2/3J), (34) 

±J 

where X 0 is the modified Bessel function. 

The physical interpretation of these equations is sim¬ 
ple. At equilibrium, the system has a finite density 
of quasiparticle n^ qp \k) and quasihole n^ qh \k) excita¬ 
tions, depending on the momentum k and simply given 
by Eqs. (30,31) which give in the small J/U limit: 

n (qp)(^) = e -P(U-4J cosk-n) 

n^ h \k) = e~^~ 2Jcosk \ (35) 

so that one simply has: 

= 1 + f [n^ qp \k) — n^ qh \k)] dk/2ir. (36) 

L J-TT 


The chemical potential /i that ensures unit filling is 
thus obtained from the (N) = L constraint: 


_ U , 1 1 2o(2/3J) 

M 2 + 2/3 n Jo (4 W 


(37) 


Similarly, one can easily compute the variance of the 
local occupation number, relying on Eqs. (25),(30),(31): 
At lowest order in J/U, we obtain: 


= 1 + 8^ + 3e-^ u ~^l 0 (4/3J) - e-^J 0 (2/3 J), 
L U z 

(38) 

which, for the value of /i at unit average filling, gives: 

= 1 + sT + 2e- /3c/ / 2 v /^o(4/3J)X 0 (2/3J). (39) 


In the limit of zero-temperature, one recovers Eq. (27) 
with only quantum fluctuations due to J. In the limit of 
vanishing J, the modified Bessel functions tend to unity, 
and one has purely thermal fluctuations. In the pres¬ 
ence of both quantum and thermal fluctuations, the vari¬ 
ance is not perfectly additive, because of crossed thermal- 
quantum effects. 


Var (tii) = (nf) - (n*) 2 

T 2 / 4 + 7r AU 

= 8 fp+J ^ [n^\k) + n^(k)} —.(41) 

which shows that quasiparticles and quasiholes equally 
contribute to the increase of the variance. 

Other expectation values can be computed as well. For 
example, the energy density is the expectation value of 
the operator Hi , Eq. (1), and is given by: 

{H i ) = -1- + Uj_ n (qp) (fc) —. (42) 

It must be emphasized that all these results are 
valid only when the density of excitations is low, 
n( qp )(£;),n( qh )(&) <C 1, i.e. when both (r^) and (n?) are 
close to unity. 

We have checked using DMRG and Quantum Monte 
Carlo numerical simulations with ALPS [26] the validity 
of Eq. (39), both in the grand canonical ensemble and in 
the canonical ensemble for large systems. For example, 
for J = 0.05U, the prediction is that the variance of the 
occupation number at zero temperature should be 0.02 
per site, while the numerical result is 0.0198453. The sit¬ 
uation is a bit more complicated for thermal excitations, 
for two independent reasons: 

• Although J/U = 0.05 seems a rather small value, 
the numerical prefactors in e.g. Eq. (20) are such 
that one has to compare 3 J with U/ 2. For J/U > 
3/2 — y/2 ~ 0.086, the Bogoliubov eigenmodes 
become unstable and the whole perturbative ap¬ 
proach breaks down. Thus, even at J/U = 0.05, 
significant corrections are expected for the thermal 
fluctuations. For f3U = 10, we find the variance 
of the occupation number to be 0.0429 while the 
prediction is 0.0385, a 10% difference. At smaller 
J/U, the relative difference is smaller. 

• Finite size effects. We observed that the variance 
obtained using the canonical ensemble in a system 
of size L has a rather strong dependence on U, 
much larger than for the grand canonical ensem¬ 
ble, at least for small J and low temperature. We 
believe that, when working in the canonical ensem¬ 
ble, fluctuations in the local occupation number are 
due to particles or holes jumping from other sites. 
If the number of sites is very large, they can pro¬ 
vide at no cost the additional particles or holes. 
In contrast, if the number of sites is too small to 
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provide these additional particles, the fluctuations 
and thus the variance will be reduced. As a rule 
of thumb, this effect is important when the prod¬ 
uct of the number of sites L by the variance of the 
occupation number is not much larger than unity. 
For L = 100, J/U = 0.05, /3U = 10, the variance is 
reduced to 0.0327. 

In order to take into account these effects, a simple 
method is to slightly reduce the density of the quasipar¬ 
ticle and quasihole excitations, Eq. (35), by multiplying 
them by a constant factor (independent of fc) slightly 
smaller than unity, which reduces the thermal fluctua¬ 
tions, leaving the quantum fluctuations unaffected. This 
reduction factor is chosen to reproduce the initial num¬ 
ber variance (at t = 0) in each sub-sample: for the ’’sharp 
gluing” scenario where the right and left sub-samples are 
initially not connected (see e.g. Fig. 8), two different 
reduction factors are used. For the ’’smooth gluing” sce¬ 
nario with an initial inhomogeneous temperature profile 
(see Fig. 10), thermal fluctuations can be globally pro¬ 
vided by both sub-samples and a single reduction factor 
is consequently used. 


E. Transport in non-equilibrium systems 

We now consider the situation of systems, not at the 
thermal equilibrium, where we want to compute the tem¬ 
poral evolution of expectation values of local operators as 
well as the currents flowing inside the system. This is a 
very complicated problem in general; we will restrict to 
the Bose-Hubbard model within the simplified assump¬ 
tion discussed above, namely when only local occupation 
numbers 0, 1 and 2 are allowed on each site. If we further 
assume that the system is everywhere close to the Mott 
insulator state with unit filling, so that the density of ex¬ 
citations is small, we can use the previous Bogoliubov de¬ 
scription in terms of quasiparticles and quasiholes. Note 
that this does not require the system to be locally at a 
thermal equilibrium. 

The Bogoliubov theory described in the previous sec¬ 
tions explicitly uses the translational invariance in order 
to obtain uncoupled elementary excitations with a well 
defined momentum. This is well suited to describe how 
excitations propagate, as it boils down to the dispersion 
relation of the excitations, Eqs. (15,17). For a system 
where translational invariance is broken by say a tem¬ 
perature gradient, this is less convenient, as one needs 
to build ’’wavepackets” coherently superimposing excita¬ 
tions with various k. A mixed position-momentum rep¬ 
resentation is in such a case more convenient, similar to 
the Wigner phase-space representation used to describe 
an ordinary quantum particle obeying the Schrodinger 
equation [36]. The temporal evolution of the Wigner 
representation is well approximated by the classical dy¬ 
namics in the semiclassical limit where the wavelength is 
much shorter than the typical size over which the Wigner 


function varies. In our case, the typical wavelength of Bo- 
goliubov excitations is 2it/k, of the order of the lattice 
spacing. Thus, provided we consider a spatially smooth 
profile of temperature/excitations, we can use a classical 
description of the system in terms of densities of quasi¬ 
particles/quasiholes n( qp )/n( qh ) depending on both the 
position x (a continuous variable in this approximation) 
and the momentum fc, which obey the classical Liouville 
equations of evolution under the Hamiltonian H e ffi 


dn( qp ) (x, fc, t) dE^^ dn( qp \x, fc, t) 
dt dfc dx 

dn( qh \x : k,t) dE)^ dn( qh \x, fc, t) 

dt = dfc dx ^ ' 

As various k components do not interact, the evolution 
is straightforward: 

7i( qp ) (X, fc, t ) = 7l( qp ) (x — fc, 0^ 

n (qh) (x, k, t ) = n (qh) (x - v { * h) t, k, o) (44) 

where the velocities of quasiparticles and quasiholes are 
given by: 


,(qp) 


dE\ 


(qp) 


d k 


(qh) 

vr = 


d E 


(qh) 


dfc 


(45) 


In the limit of small J, they are, see Eqs. (15) and (17): 


= 4 J sin k 


„( qh ) 


2 J sinfc, 


(46) 


which means that quasiparticles propagate twice faster 
than quasiholes. 

A key and non-trivial point of this approach is that it 
predicts that quasiparticle/quasiholes excitations propa¬ 
gate ballistically , as suggested by the quasi-exact numer¬ 
ical results using METTS presented above. 

Eqs. (44) make it trivial to compute the total density 
of quasiparticles and quasiholes at a given position: 



n (qh) (*, t) = [ +W ^ n (qh) (x - 4 qh k k, 0) (47) 

J-7T 27T V J 


From these formula, one can deduce the expectation val¬ 
ues of any local observable, following the derivation in 
section HID: 


(n(x)} = 1 + n( qp ) (x, t) — n( qh ) (x, t) 

(n 2 (x)) = 1 + 8 ^ + 3 n^ p \x, t ) — n^ qh \x, t) 

Var (n(x)) = 8^- + n^ qp ^ (x, t) + n^ qh ^ (x, t) 

(ff(x))/C/=-4T +n (qp)( X;i ) ( 48 ) 
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Although the quasiparticle/quasiholes densities are not 
in principle directly observable, these equations show 
that it is enough to measure two independent quantities, 
for example the average local occupation number and its 
variance to extract both the quasiparticle and the quasi¬ 
hole densities. 

The Bogoliubov approach also makes it possible to 
compute the currents. From Eq. (43), it follows that 
one can define quasiparticle/quasihole currents: 

j( qp )(x, fc, t) = n^ qp \x,k,t) 

j( qh )(x, &, t) = n( qh \x,k,t) (49) 

which automatically satisfy the continuity equation 
dn/dt + dj/dx = 0. 

The current at a given position is simply summed over 
all momenta contributions: 

j i( qp) (x,t) = [ v i qp ^ n^ qp ) (x, fc, t) — 

J — 7T 27r 

i (qh) (*, t) = J+" 4 qh) (x, k, t ) g (50) 

Following Eqs. (48), the mass and energy currents are 
thus given by: 

j n (x,t)=j^\x,t)-j^ h \x,t) 
j E (x,t)/U=j^\x,t). (51) 

The physical interpretation of these equations is clear: in 
the small J limit, the dominant contribution to energy is 
the interaction brought by sites with double occupation, 
associated with quasiparticle excitations. In contrast, a 
quasihole does not lead to a change of local energy and 
thus does not contribute to the energy current. As ex¬ 
pected, quasiparticles and quasiholes both contribute to 
the mass (density) current, with opposite signs. 

The previous set of equations describe the ballistic 
transport of quasiparticles and quasiholes in the system. 
In general, the momentum distribution of quasiparticles 
and quasiholes at a given position is not given by a ther¬ 
mal distribution, Eq. (35). This implies that our descrip¬ 
tion goes beyond a local thermal equilibrium. Even if 
the initial state at t = 0 is in a local thermal equilibrium 
described by Eq. (35), for space dependent (inverse) tem¬ 
perature /3(x) and chemical potential /i(x), this property 
is lost during time evolution. 

While the integrals involved in Eqs. (47) have a trivial 
structure and are easily numerically computed for arbi¬ 
trary initial distributions, it is in general difficult to per¬ 
form the integrals analytically. There are however simple 
cases where it is possible. 

Let us first consider the situation where two half 
blocks, each at thermal equilibrium with unit filling, but 
with different temperatures (/ $l,Pr ) (and consequently 
different chemical potentials given by Eq. (37)) are con¬ 
nected at time t — 0. Taking x = 0 as the connection 


point, this implies that the initial quasiparticle/quasihole 
distributions are given by: 

n( qp ) (x, &, 0) =q~Pl(U- 4J cosk-fi L ) f or x < 

n( qp ) (x, k, 0) = e -P R ( u - qj cosk-n R ) f or x > 

n^ qh \x : k, 0) = cosk) for x < 0, 

n (qh) (x,/c,0) = e -P R ^ R - 2Jcosk ) for x > 0. (52) 


The temporal evolution of each k component is trivial: 
it consists of a step function moving at a velocity 

(resp. ?4 qh )) f° r quasiparticles (resp. quasiholes). In 
effect the densities have a very simple geometrical struc¬ 
ture: they are ’’smoothed” steps connecting the asymp¬ 
totic ’’left” density (on the left side) to the asymptotic 
’’right” density (on the right side). This step function is 
centered around the origin, keeps the same shape at any 
time, being simply stretched along the x-axis proportion¬ 
ally to time. Because quasiparticles (resp. quasiholes) 
have a maximum velocity 4J (resp. 2J), the smooth 
step extends only inside a ’’light cone”, in a finite range 
[—4Jt, 4Jt] (resp. [—2Jt,2Jt]). The situation is simi¬ 
lar for the currents which vanish outside the light cones 
and keep the same spatial structure, simply stretching 
linearly in time. 

We could not perform analytically the integral over k 
for the density, but could evaluate it for the currents: 


j (qp \x,t) = - 
7r 

j^\x,t) = - 

7 T 


Pl Pr 


'-Plu. Sinh /^ _ Sinh ^ 

Pl Pr 


(53) 


where u = yi6./' 2 — x 2 /t' 2 . v = y4./ 2 — x 2 /t 2 . These ex- 
pressions are valid inside the ’’light cones” \x\ < 4 Jt (resp. 
\x\ < 2 Jt) for the quasiparticles (resp. quasiholes); out¬ 
side the light cones, the currents vanish. 

The currents at the origin x = 0 are especially simple. 
They do not depend on time and are given: 


j (qp) ( 0) = - 

7T 

i (qh) (o) = - 

7 T 


P -faT/-^) sinh4 ^ J _ p -/w-M B ) sinh4 ^ J 

Pl Pr 


'-PM sinh2 ^ J _ p -to smh2 p R J ' 
Pl Pr 


(54) 


The fact that they are time-independent is a signature 
of ballistic propagation of the excitations, and nicely fits 
the observations made in Fig. 4 on our quasi-exact nu¬ 
merical METTS simulations. Both currents are simply 
of the form #(/3 l) — #(/Jr) for some function g , a property 
already emphasized as a signature of ballistic transport 
in [15]. In contrast, diffusive transport would be charac¬ 
terized by currents depending on the local gradients. 

In the regime of intermediate temperature J 
^ < I/, the expression simplifies and the currents 
have a ’’semi-circle” spatial shape: 


j (qp) (x,t) OC y/l6J 2 — X 2 ft 2 

j( qh \x,t) OC y/4J 2 — X 2 /t 2 . (55) 
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At lower temperature (/3J of the order of unity or 
smaller), the shape is qualitatively similar. 


IV. COMPARISON OF THE BOGOLIUBOV 
THEORY WITH NUMERICAL DATA 

A. Sharp gluing 

We now give a more detailed analysis of the simulations 
already presented in Fig. 1 and Figs. 3-5. Recall that we 
consider two Mott insulators of length M — 35, each one 
at thermal equilibrium, but with different temperatures 
/3lU = 6, [3rU = 8. Each insulator is represented by an 
appropriate thermal density matrix (METTS ensemble), 
the total state of the system being the tensor product 
of those density operators. At time t = 0 the tunneling 
between touching sites is turned on to a common value 
J/U = 0.05. While in Figs. 1 and 3 we have shown the av¬ 
erage particle density (n^), squared particle density (nf) 
and variance Var(n^) = (n?) — (n*) 2 at selected times, 
Fig. 6 presents the same data resulting from a quasi- 
exact numerical simulation as color plots (left column) 
together with the predictions given by the Bogoliubov 
theory (right column) described in the previous Section. 

A good agreement between the quasi-exact METTS 
simulations and the Bogoliubov prediction is visible al¬ 
ready in the top row, where we compare densities, (n^): 
a bump propagating to the right side, a hole propagating 
to the left side. On the left side, the higher temperature 
implies a higher density of quasiparticles and quasiholes. 
Once the two parts are connected, the quasiparticles 
propagate ballistically at twice the speed of the quasi¬ 
holes. It results in an excess of quasiparticles - hence a 
higher density - between the quasihole light cone (dashed 
black line) and the quasiparticle light cone (white dashed 
line). As noted earlier, the variance Var(r^), shown in 
the second row, is less sensitive to noise, and a very good 
agreement between numerical METTS simulations and 
the Bogoliubov theory is observed. 

While the basic ingredients of the theoretical approach 
are the quasiparticle and quasiholes densities, these quan¬ 
tities cannot be directly measured in the METTS simu¬ 
lations. It is, however, possible to invert Eqs. (48) and to 
deduce the densities from the measured quantities ( rii ) 
and (nf): 

(qp) _ Var(n^) -1- ( tl ^) 1 4 J 

~ 2 W 

(qh) _ Var(rij) — (n*) + 1 4 J 2 

n i ~ 2 W ( j 

These formulae allow us to determine the distribution 
in time of quasiparticles and quasiholes for simulations, 
as well as conversely find the predictions for the particle 
density and its variance from the quasiparticle distribu¬ 
tions. Clearly the excess of particles is observed on the 
cold side (with the corresponding hole on the hot side due 


to conservation of the total number of particles) moving 
out ballistically from the junction point. That is due to 
the spread of quasiparticles, as described in the previous 
section. 

Observe that while we keep the same scale for the den¬ 
sity, we use different scales for the variance (second row in 
Fig. 6) and for the quasiparticle (third row) and quasihole 
(fourth row) densities. That is due to finite size effects 
discussed extensively in the previous section. Instead of 
fitting the correction we have chosen to plot the raw data; 
the correction is “automatically” taken into account by 
the rescaling of the colorbar. This provides a convincing 
argument that a single correction factor valid at all times 
is sufficient to bring the results of METTS simulations 
and the theoretical predictions together. While of course 
some differences may be visible, the overall agreement is 
quite remarkable showing beyond doubt that the Bogoli¬ 
ubov theory describes very well the results of the simu¬ 
lations and thus the physics of heat and mass transfer in 
this system. Note especially the clear demonstration that 
quasiholes propagate more slowly than quasiparticles (by 
a factor 2), as the quasihole density is affected only inside 
the inner light cone - marked with the black dashed lines 
- while the quasiparticle density is affected inside the 
outer light cone - marked with the white dashed lines. 

These conclusions are further confirmed by inspection 
of different possible currents. As before, the simulations 
give us access to mass current j n , Eq. (8) as well as to the 
energy current j E , Eq. (6). Those can be related to quasi¬ 
particle and quasihole currents via Eqs. (51). One can 
also define the variance current, which following Eq. (48) 
can be defined as 




j n (x,t) + 


2j g QM) 

U 


(57) 

Since the variance is directly related to temperature, the 
variance current can be considered as a heat current. 

Observe in Fig. 7 that, as expected, the currents are 
nonzero only in the “light-cone” region around the point 
of junction. The agreement between the simulations and 
the theory is less spectacular for mass and energy cur¬ 
rents while it is much better for the heat current and, in 
particular, quasihole current (bottom row). Indeed the 
cone both for simulations and theory is then restricted 
by the maximal velocity of quasiholes. The differences in 
results between the two approaches may be again traced 
back to finite size effects. 

Color plots in Figs. 6,7 show that our Bogoliubov ap¬ 
proach is qualitatively correct: heat and mass transport 
in the system are conveniently described by quasiparti¬ 
cle/quasihole excitations that propagate ballistically. In 
order to make the comparison more quantitative, we show 
in Fig. 8 the average particle density (rii) and the average 
squared density (n?) at time tU = 60. After the proper 
rescaling necessary due to finite size effects (see discus¬ 
sion above), the agreement is very good. The density 
bump (on the right side) and hole (on the left side) is 
well predicted. For the average squared density (nf), 
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FIG. 6. (color online) Comparison of quasi-exact numerical METTS simulations (left column) with the Bogoliubov theory 
(right). From top to bottom: mean site occupation (ni), variance (n 2 ) — (n*) 2 , quasiparticle density n[ qp ^ and quasihole density 
n[ qh \ The dashed white (resp. black) straight lines correspond to ballistic propagation of quasiparticles (resp. quasiholes) from 
the connection point between the two subsamples at the maximal allowed velocity, i.e. x^ qp \t) = ±4 Jt (resp. x^ qh \t) = ±2Jt). 
See text for a detailed discussion. 
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FIG. 7. (color online) Comparison of the currents computed from the numerical METTS simulations (left column) with the 
Bogoliubov theory (right). The currents are obtained by numerical calculations for the same parameter values as in Fig. 6. 
From top to bottom: mass current, energy current, heat (variance) current as defined in Eq. (57) and quasihole current as 
defined in Eqs. (51) (the quasiparticle current coincides with the energy current). See text for discussion. 
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FIG. 8. (color online) Comparison of the average particle 
density (m) (lower curves) and the average squared density 
(n?) (upper curves) at time tU = 60, for the quasi-exact 
METTS simulations (solid black lines) and the Bogoliubov 
theory. The raw Bogoliubov predictions are the dotted green 
curves; rescaled results, taking into account finite size effects, 
as shown as dashed red lines and agree very well with the 
METTS simulations. See text for discussion of the various 
features. 
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FIG. 9. (color online) Comparison of the mass and energy 
currents for the same data than in Fig. 8. The energy cur¬ 
rent is predicted by the Bogoliubov theory to depend only 
on the quasiparticle current, see Eq. (51). It displays the 
characteristic ’’semi-circle” shape predicted by our Bogoliubov 
theory, Eq. (55). In contrast, the mass current is the differ¬ 
ence between the quasiparticle and quasihole ’’semi-circles”, 
see Eq. (51), and presents two maxima propagating at the 
maximum velocity of the quasiholes ±2 J. 


one can see three distinct regions between the ’’hot” left 
plateau (not yet affected) and the ’’cold” right plateau: 
two rather abrupt cliffs on the edges of the two plateaus 
separated by an intermediate much flatter region. The 
two cliffs correspond to regions already reached by quasi¬ 
particles while the intermediate flat region is affected 
both by quasiparticles and quasiholes. The kinks at 
the frontiers between regions correspond to quasiparti¬ 
cles/quasiholes with maximum velocities 4J/2J, respec¬ 
tively. They are smoothed out, but still clearly visible, 


in the quasi-exact METTS simulations. Note that all 
quasiparticles/quasiholes in the system are thermally ex¬ 
cited, so that the effects observed are truly due to non¬ 
equilibrium thermodynamical properties. 

In Fig. 9, we show the comparison between the mass 
and energy currents as computed from the METTS sim¬ 
ulations and the predictions of the Bogoliubov approach. 
Not surprisingly, the salient features are quantitatively 
well predicted: the energy current, predicted by Eq. (51) 
to depend only on the quasiparticle current displays a 
single bump with the characteristic ’’semi-circle” shape 
predicted by Eq. (55); the mass current is the differ¬ 
ence between the quasiparticle and quasihole currents 
and consequently reveals two maxima propagating at the 
maximum velocity of quasiholes ±2 J. 


B. Smooth gluing 

We now confront the predictions of our Bogoliubov 
theory with the numerical data obtained for a smoothly 
glued sample, as the one described in section IIB 2. As 
previously, the whole system consists of 70 sites with 
PlU = 6 on the hot side and PrU = 8 on the cold 
side. Fig. 10 shows the variance of the particle density 
across the system at initial time as well as after some 
real time evolution. In order to compare the theory and 
the numerical experiment correctly, one has to introduce 
some rescaling, taking care of the finite size effect. While 
two different rescaling factors were used for the ’’sharp 
gluing” case where the two subsamples are initially not 
connected, we here use a single rescaling parameter for 
the whole sample, and all times. This is because - by con¬ 
struction - a smooth gluing procedure find the equilib¬ 
rium initial state for the whole sample at the initial time. 
As before, real time evolution washes out site to site fluc¬ 
tuations of the variance: starting from, say tU = 10, the 
variance is perfectly reproduced by the Bogoliubov the¬ 
ory. 

Fig. 11 presents the mass and the energy currents at 
tU = 60 using the same single scaling factor. Clearly 
some quantitative discrepancies between the numerics 
and theoretical simulations exist. One may, however, ob¬ 
serve quite nice qualitative shape agreement between the 
two curves. 

As a final example, we consider the same system but at 
lower temperature. This time the hot part corresponds 
to PlU = 8 while the cold part to PrU = 10, again with 
a smooth gluing of the two subsystems. The results are 
presented in Fig. 12 and Fig. 13. While the fluctuations 
are notably larger (quantum effects are more visible), the 
behavior of the system is similar to the previously con¬ 
sidered case. Note also that the imperfect quantitative 
agreement for the currents observed at higher temper¬ 
ature, Fig. 11, is here significantly better. A possible 
explanation is that, at higher temperature, our simple 
Bogoliubov theory breaks down at shorter times, either 
because of additional excitations or because of the finite 
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FIG. 10. (color online) Particle number variance at initial 
time t — 0 and at tU = 60, for two smoothly glued Mott in¬ 
sulators with (inverse) temperatures /3lU = 6, /3rU = 8, and 
J/U — 0.05. The interface size, Eq. (2), is s = 15. Solid lines 
correspond to quasi-exact numerical results using the METTS 
method, while dashed curves are theoretical estimates. A sin¬ 
gle parameter - required to take into account finite size ef¬ 
fects - has been used to adjust the numerical and theoretical 
curves. The sharpest (black) step corresponds to the initial 
time t = 0. The smoother (red) curve corresponds to tU = 60. 
Note that the real time evolution strongly smooths the initial 
site-to-site fluctuations. At tU = 60, the theoretical predic¬ 
tion is almost indistinguishable from the quasi-exact numeri¬ 
cal simulation. 



FIG. 11. (color online) Mass (upper panel) and energy (lower 
panel) currents at tU = 60 for the same data as in Fig. 10. Nu¬ 
merical results (solid black curves) are compared with the pre¬ 
dictions of the Bogoliubov theory (dashed red curves). While 
the magnitude of the currents is not predicted properly, the 
overall shape is faithfully predicted by the theory. 

lifetime of quasiparticle/quasihole excitations. 


V. CONCLUSIONS 

We have studied non-equilibrium dynamics occurring 
when two quantum insulators at different temperatures 
are brought into contact. The process, in the long run, 
should lead to equilibration of temperatures (at least in 
the thermodynamic limit). Using Mott insulators of the 



FIG. 12. (color online) Particle number variance at initial 
time t = 0, tU =10 and at tU = 60, for two smoothly glued 
Mott insulators with (inverse) temperatures /3lU = 8 ,/3rU = 
10, and J/U = 0.05. The interface size, Eq. (2), is s = 15. 
Solid lines correspond to quasi-exact numerical results using 
the METTS method, while dashed curves are theoretical es¬ 
timates. A single parameter - required to take into account 
finite size effects - has been used to adjust the numerical and 
theoretical curves. The sharpest (black) step corresponds to 
the initial time t — 0. The smoother (red and green) curves 
corresponds to tU = 10 and tU = 60. The real time evolu¬ 
tion strongly smooths the initial site-to-site fluctuations. At 
tU = 60, the theoretical prediction is almost indistinguishable 
from the quasi-exact numerical simulation. 



FIG. 13. (color online) Mass (upper panel) and energy (lower 
panel) currents at tU = 60 for the same data as in Fig. 12. 
Numerical results (solid black curves) are compared with the 
predictions of the Bogoliubov theory (dashed red curves). The 
agreement is very satisfactory. 

Bose-Hubbard model as a specific example - as it is close 
to experimental realization with ultracold atoms - we, 
however, observe that heat transport is initially ballistic: 
it efficiently transfers energy from the hot to the cold side 
of the system, but the full thermalization may occur only 
on a longer time scale, not reachable by our numerical 
simulations. 

Fully numerical, “quasi-exact” calculations have been 
performed using the METTS algorithm [18]. Even with 
important computer resources, we could only study the 
dynamics on a rather limited amount of time. This is 













16 


due to two primary reasons. First, the METTS ensemble 
needed to obtain reasonable averages consisted of more 
than 10 000 vectors; each of them had to be evolved us¬ 
ing the TEBD (or t-DMRG) algorithm. The price to 
pay is a direct 10 4 increase of the computational time in 
comparison to T = 0 studies of e.g. quantum quenches. 
Second, for T > 0, the METTS ensemble includes by con¬ 
struction some significantly excited states; the real time 
evolution of such states leads to a fast increase of the 
entanglement, requiring significant extension of the bond 
dimension in the MPS description. A posteriori, one may 
conclude that this is not the method of choice and other 
approaches such as a purification scheme [15] or meth¬ 
ods employing Heisenberg picture evolution [IT] could 
perform much better, at least according to comments in 
these papers. A very recent comparison of purification 
schemes with METTS [22] reached similar conclusions. 

The numerical results have been quantitatively com¬ 
pared with a simple analytic theory based on Bogoliubov 
excitations. The theory explains the ballistic transport 
numerically observed in our system providing good esti¬ 
mates for observable quantities, provided some plausible 
corrections are included that took into account the finite 
size effects. The theory applies also for smoothly glued 
system: the transport may be fully understood as the 


movement of quasiparticles and quasiholes of the inte¬ 
grate Bogoliubov theory. In some sense, this explains 
the success of the METTS approach for description of 
this particular system. As conjectured by Prosen [16], 
the entanglement growth in an integrable system is slower 
and easier to handle. 
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